This notebook presents an agent-based model in combination with a SIR-like compartmental model of powdery mildew infection spread among grapevines in a spatial grid.
This is a study extending on a simple exploration done by Paul Petersik, with added details on consideration of infection and recovery transition and different aspects reflecting spatial heterogeneity displayed in the process of fungal infection spread.
Modelling and visualisation is done in python with AgentPy.
# Model design
import agentpy as ap
# Visualization
import matplotlib.pyplot as plt
import IPython
# General Utility
import random
This section details the models used in this study and a guide to walk through the python implementation of the model.
This study divides the grapevine population into 4 compartment:
Susceptible
Latent
Infected
Recovered
A class AgentState is created to host the 4 states/compartments for the model.
class AgentState():
SUSCEPTIBLE = 0
LATENT = 1
INFECTED = 2
RECOVERED = 3
The SLIR (Susceptible-Latent-Infected-Recovered) model runs inside a cellular automaton, with a collection of agents spatially represented as a grid. In context, the grid represents a sample field of grape vines and each cell, hosting an agent, represents an individual grape vine. As the automaton runs and iterates over generations, the agents updates their status according the pass of time and their surrounding circumstances.
The automaton has two main tasks: initialisation and stepping through generations.
The setup function is run when the model/automaton gets initialised. It prepares the grid and its agents, before the model starts to run, with an initial generation. It starts by creating the specified number of agents in this model and assigns them to the grid, filling on possible cells. The agents are then prepared with corresponding initial attributes, differing according to their spacial distribution. Such differences reflects spacial heterogeneity in grape vines' growth condition and inert ability.
The grapes vines are initialised as either part of the susceptible population or the resistant population. Such difference will be explained in the next section with the logic where spacial heterogeneity is introduced. However, both these compartments represent healthy grape vines. In order to set the automaton running, some infected individuals are randomly planted into the grid. These infected vines are in latent state so that their symptom will develop gradually until active and start to infect neighbouring vines, spreading the disease.
The model moves forward in time through generations of agents by calling the step function. The method is called in every iteration of generation, reflecting changes happening to grapevines in the vineyard. The function determines the next state of each individual agents in the grid, by first finding its neighbourhood. The state of a grapevine is influenced by infectious state of its surrounding vines in the spreading process of powdery mildew infection. The specific state-transition flow is detailed in the Agent section.
Once the next state of each agent is determined, the model sets the next state of each agent as its current state, signifying that a generation has passed and the necessary state transitions happening during this period have been done.
The model stops before reaching the maximum number of iterations when there are no more latent or infected agents in the vineyard. Since without these agents as active players, all agents in the grid remains static and simulation virtually ceases, no matter how time passes.
model_generation is a global variable that is updated by the automaton and tracks the current timestamp of the simulation. This variable is used later for time-related agent state transitions.
model_generation = 0
class VineInfectionModel(ap.Model):
def setup(self):
# Create agents
vines = self.agents = ap.AgentList(self, self.p.size**2)
# Create grid and add agents to the grid
self.vineyard = ap.Grid(self, (self.p.size, self.p.size))
self.vineyard.add_agents(vines)
# Initiate a state variable for all agents
set_up_agent_variable_attributes(self.vineyard)
# Initialise model with infected vines
initial_infected = self.vineyard.agents.to_list().random(
self.p.initial_infection_count)
initial_infected.state = AgentState.LATENT
def step(self):
global model_generation
# Calculate the next state of each agent
for vine in self.agents:
neighboring_vines = self.vineyard.neighbors(vine).to_list()
update_agent(vine, neighboring_vines)
# Update the state of all agents
for vine in self.agents:
vine.state = vine.state_next
# Update model generation count
model_generation = self.t
# Stop simulation if no active agent is left
latent_agent_count = len(self.agents.select(self.agents.state == AgentState.LATENT))
infected_agent_count = len(self.agents.select(self.agents.state == AgentState.INFECTED))
if (latent_agent_count + infected_agent_count) == 0:
self.stop()
Spatial heterogeneity is a main area of focus of this study. Measures are introduced to this simulation to simulate and reflect how the spread patterns of infection is affected spatially due to various factors.
In the current implementation, some of the main influencer of spread are:
The entire grid of vinegrapes are intercropped with a cultivar with more resistance to catch infection. This is done in order to setup barriers in between groups of normal grapevines to break the wave of fungal spread and slow down disease propagation.
The resistant cultivar perform much better against the infectant. However, it is not realistic to have too many such individuals in the overall population, as mentioned by Paul Petersik.
So, you may now ask the question: “Why not only plant resistant varieties?” The problem is, most consumers do not really like the wine made from resistant varieties (yet).
The following function demonstrates how intercropping is implemented in this model. The function uses the coordinates of the agent to determine whether the agent of interest is a normal grapevine cultivar or an intercropped resistant cultivar. Note that the pattern of intercropping is a evenly distributed row-by-row setup. As a result, whether an agent is intercropped resistant cultivar can be deduced by taking the remainder of the total pattern distance.
The above image visualises the intercropping setup, where light green spaces represent the normal grapevine cultivar and the darker green spaces represent the resistant cultivar.
intercropping_distance = 5
intercropping_row_count = 5
def is_intercropped_resistant_cultivar(position):
x_coordinate = position[0]
y_coordinate = position[1]
return x_coordinate % (intercropping_distance + intercropping_row_count) >= intercropping_distance
default_infectivity_offset = 10
resistant_infectivity_offset = 180
default_latency_duration = 8
def set_up_agent_variable_attributes(grid: ap.Grid):
for agent in grid.agents:
position = grid.positions[agent]
is_resistant_cultivar = is_intercropped_resistant_cultivar(position)
# Set cell state
agent.initial_state = AgentState.RECOVERED if is_resistant_cultivar else AgentState.SUSCEPTIBLE
agent.state = agent.initial_state
# Set infectivity offset based on cultivar
agent.infectivity_offset = resistant_infectivity_offset if is_resistant_cultivar else default_infectivity_offset
# Set latency duration
agent.latency_counter = 0
agent.latency_duration = random.gauss(default_latency_duration, 0.5)
In the cellular automaton, each agent of the model functions as a state-machine that updates or maintains its state while the automaton's generation advances. As explained in the previous section, this is achieved by progressively calling the update_agent function to calculate the next state of every agent in the grid, based on its current state and its surrounding neighbors.
Putting the state machine in context, each vine (agent) can belong to a compartment shown in the previous section, and can transition into another compartment as it get infected by powdery mildew.
def update_agent(vine: ap.Agent, neighbors: ap.AgentList):
if vine.state == AgentState.SUSCEPTIBLE or vine.state == AgentState.RECOVERED:
vine.state_next = AgentState.LATENT if s_2_l(vine, neighbors) else vine.state
if vine.state == AgentState.LATENT:
vine.latency_counter += 1
vine.state_next = AgentState.INFECTED if l_2_i(vine, neighbors) else vine.state
if vine.state == AgentState.INFECTED:
if self_recovery(vine, neighbors):
vine.initial_state = vine.state_next = AgentState.RECOVERED
vine.infectivity_offset = resistant_infectivity_offset
return
if manual_removal(vine, neighbors):
vine.state_next = vine.initial_state
return
vine.state_next = vine.state
This section describes the process of each state transition of the aforementioned state machine. The
default_infectivity = 10
def s_2_l(vine: ap.Agent, neighbors: ap.AgentList):
neighborhood_infectivity = default_infectivity * len(neighbors.select(neighbors.state == AgentState.INFECTED))
prob = random.random()
return prob <= (neighborhood_infectivity / (neighborhood_infectivity + vine.infectivity_offset))
def l_2_i(vine: ap.Agent, neighbors: ap.AgentIter):
return vine.latency_counter >= vine.latency_duration
removal_threshold = 10
recovery_rate = 0.001
# TODO: make more sense out of this
def self_recovery(vine: ap.Agent, neighbors: ap.AgentIter):
return random.random() <= recovery_rate
def manual_removal(vine: ap.Agent, neighbors: ap.AgentIter):
return model_generation > 0 and model_generation % removal_threshold == 0
parameters = {
'size': 100,
'initial_infection_count': 2,
'steps': 200,
}
# Create single-run animation with custom colors
def animation_plot(model:VineInfectionModel, ax):
attr_grid = model.vineyard.attr_grid('state')
color_dict = {
0:'#7FC97F', # Susceptible
1:'#e6bc2c', # Latent
2:'#d62c2c', # Infected
3:'#1FC97F', # Reovered / Resistant
None:'#d5e5d5'
}
ap.gridplot(attr_grid, ax=ax, color_dict=color_dict, convert=True)
ax.set_title(f"Simulation of a SIR model\n"
f"Time-step: {model.t}, infected vines: "
f"{len(model.agents.select(model.agents.state == AgentState.INFECTED))}")
fig, ax = plt.subplots()
model = VineInfectionModel(parameters)
animation = ap.animate(model, fig, ax, animation_plot)
IPython.display.HTML(animation.to_jshtml(fps=15))